Path Integral Calculations of the Hydrogen Hugoniot Using Augmented Nodes 
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We calculate the hydrogen Hugoniot using ab initio path integral Monte Carlo. We introduce 
an efficient finite-temperature fixed-node approximation for handling fermions, which includes an 
optimized mixture of free particle states and atomic orbitals. The calculated Hugoniot confirms 
previous fixed-node path integral calculations at temperatures around T = 30 000 K and above, 
while approaching smoothly the low temperature gas gun results. The ability to optimize the free 
energy within the path integral opens many new possibilities for developing nodal density matrices 
for path integral simulations of other chemical systems. 
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Predicting the properties of materials at extreme con- 
ditions of high pressures and high temperatures is a fun- 
damental challenge in computational physics, particu- 
larly when quantum effects are significant. One such sys- 
tem is the hydrogen Hugoniot, measured in laser shock- 
wave [TJ [5] and gas gun [3] experiments, which has be- 
come one of the benchmarks for state-of-the-art theory. 

Path Integral Monte Carlo (PIMC) allows direct calcu- 
lation of thermodynamic properties of many-body quan- 
tum systems at finite temperature [4HT0]. For fermions, 
an additional "fixed-node" constraint is imposed to pre- 
vent an exponential slowdown in convergence with in- 
creasing system size and decreasing temperature [TTJ [T^] . 
This requirement to supply an a priori form for the nodal 
surface of the many-body wavefunction has limited the 
application of PIMC almost exclusively to hydrogen and 
helium systems at high temperatures. Not only does 
the use of an approximate nodal function introduce un- 
controlled error, traditional algorithms for enforcing the 
nodal constraint have been shown to become increas- 
ingly inefficient at lower temperatures, ultimately failing 
to satisfy the Monte Carlo ergodicity requirement. For 
hydrogen, nodal errors begin to manifest at lower tem- 
peratures, below 100 000 K [7J, and ergodicity becomes 
difficult to achieve below 20 000 K [7J. 

In this Letter, we present a simplifying approximation 
to fixed-node PIMC that resolves several of these issues. 
This new implementation is fully symmetric in imaginary 
time, thus achieving better computational efficiency in 
addition to avoiding the ergodicity issues of the tradi- 
tional approach. More importantly, the method allows 
the relative free energies of different nodal choices to be 
sampled directly within PIMC simulations, thus provid- 
ing a systematic way to construct an optimal nodal func- 
tion. This is an important independent verification of 
other ab initio proposals and brings the benefits of PIMC 
to a much broader physical domain. While the ability 
to test and improve nodes using total energies has been 
common in ground-state Quantum Monte Carlo (QMC), 
extending practical tests of nodes to compare free ener- 
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FIG. 1. (color online) The Deuterium Hugoniot curve calcu- 
lated within PIMC from T = 125 000 K down to T = 5 000 K. 
While free particle nodes are sufficient at high temperatures 
(above 125 000 K, see Table [I]), augmented nodes are neces- 
sary once bound electronic states form. Our Hugoniot curve 
approaches the correct free particle limit of fourfold compres- 
sion in the high temperature and pressure limit (not shown). 



gies is an important advancement in the study of finite 
temperature nodes. This could have far-reaching conse- 
quences, as atomic and molecular orbitals can be intro- 
duced without bias, potentially paving the way to extend 
PIMC to heavier elements and lower temperatures. 

We demonstrate the new formalism with calculations 
of the hydrogen Hugoniot (Fig. [TJ, using augmented 
nodes incorporating Is atomic orbitals. Our calcula- 
tions of the hydrogen Hugoniot span temperatures from 
T = 1 000 000 K all the way down to 5 000 K, overlapping 
for the first time with ground state theories like Den- 
sity Functional Theory (DFT) [T31 [H] , as well as lower- 



2 




(a) reference- slice constraint (b) antipodal-slice constraint 



FIG. 2. Schematic illustration of the (a) reference-slice fixed- 
node, Eq. Q, [11] and (b) antipodal-slice fixed-node, Eq. |3|, 
constraints for path integrals. 



temperature gas gun experiments There is a general 
consensus over the high and low pressure limits of the 
Hugoniot curve. A controversy arises at the intermedi- 
ate pressure, around 1 Mbar, where laser-driven shock 
measurements pQ observe high compression 77 ~ 5.5 — 6 
while magnetically driven flyers [15] and convergent ex- 
plosives [16] find a stiffcr curve closer to 77 ~ 4. We agree 
with the latter with a maximum compression measure- 
ment t] ~ 4.2. We see no evidence of a plasma phase tran- 
sition that has been predicted by some other simulations 
[T71 ITS] . Our results are consistent with recent QMC and 
DFT calculations that show that a metal-insulator phase 
transition occurs only stakes place at lower temperatures, 
below 2 000 K, at pressures above 1.2 Mbar [lf)ll2T)] . 

PIMC uses Feynman's path integral formulation of sta- 
tistical mechanics to describe a quantum system in ther- 
mal equilibrium with a path integral over imaginary time: 

(R\e-P H \R>) = L m) = R (1) 

R{0) = R' 

where R represents all the particle coordinates, and the 
path R(r) extends for an imaginary time /3h = h/ksT. In 
PIMC simulations, the path is discretized into Nt slices, 
where Nt is the Trotter number, which divides imaginary 
time into short imaginary-time intervals, At = (3h/Nx- 
To handle fermions, the path integral, Eq. ([!]), is anti- 
symmetrized over all permutations of identical particles, 
which leads to an exponential decrease in algorithmic effi- 
ciency at low temperatures and large system sized, known 
as the fermion sign problem. 

One strategy to manage the fermion sign problem is the 
fixed-node approximation. The traditional implementa- 
tion of the fixed-node approximation restricts the path 
integral to paths that satisfy the constraint [TTj . 

Pt (R(t),R*;\t\)>0 (VrG[-^/2,W2]), (2) 

where the reference slice coordinates are taken at r = 
so that R* = i?(0), as shown in Fig.[2^a). This rigorously 
enforces the nodal constraint, but adds imaginary-time 
dependence to the evaluation of the nodal action [TT] . 
The nodal model pr must be defined for imaginary times 



ranging from Ar to /3h/2, corresponding to temperatures 
ranging from 2T to N T T. 

In the antipodal slice approach, we make a simplifying 
approximation, 

p T (R(r),R(T + Ph/2)^0 (VrG [-0h/2,O]), (3) 

illustrated in Fig. [2jb). For each slice R(t), we evaluate 
the constraint with the slice R(t+/3H/2) that is antipodal 
to the slice, rather than the reference slice R* . While we 
have tested this constraint on many problems, we have 
not determined its range of validity; the results presented 
here are a pragmatic test of the antipodal approximation. 
By lifting the imaginary-time dependence, the antipo- 
dal slice approximation allows one to avoid ergodicity 
errors caused by the reference slice getting "pinned" at 
lower temperatures where the number of imaginary-time 
slices Nt becomes large and long permutation cycles be- 
come increasingly common. The nodal model is greatly 
simplified, as pt now only needs to be evaluated at one 
temperature 2T for the imaginary time r = f3h/2. Par- 
allel scalability is similarly increased, as antipodal col- 
lections of slices can simultaneously make multiple mul- 
tilevel moves [4j on local processors without requiring 
global communication. 

The most significant advancement allowed by the an- 
tipodal approximation, however, is the ability to directly 
compare the free energies of different nodal models. We 
hypothesize that the best nodal choice will the lowest free 
energy, in analogy to the best ground state wavefunction 
nodes having the lowest energy. For the nodal density 
matrix, we take a product of Slater determinants, 

p T {R,R') =det |p(r^ t ,r it )| det |p(r 4 , r^)|, (4) 

where p(r,r) are single-particle density matrices at a 
temperature T/2. In this work we augment the single 
particle density matrices by including a hydrogen Is or- 
bital in the thermal mixture, 



P(r ' r; ^' w)= (2^/m)^ 

N lon (5) 
k 

where the weight w is a free parameter and f/ , is(r) is the 
normalized hydrogen Is orbital. Here the first term is 
the free particle thermal density matrix, and the sum 
is a mixture of density matrices for hydrogenic orbitals 
about each of the ions, each with a weight w. 

To find the optimal value of w, we histogram free en- 
ergy differences for nodes with different values of w using 
the acceptance ratio method [Hj , as shown in Fig. [3j At 
T — 62 500 K a value of w around 4.7 increases the par- 
tition function by about a factor of four. Note that the 
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FIG. 3. (color online) The relative partition function Z oc 
exp(-/3F[p T (R, R';T,w)]) for the system with r s = 2.0, for 
different nodal choices. The fixed-node free energy F[px] is a 
functional of the nodal density matrix, pr- The parameter w 
is the weight of the Is orbital in the augmented nodal mixture, 
Eq. (JsJ . At higher temperatures, there is a clear maximum in 
the free energy, while the free energy saturates for large values 
of w at low temperatures. Arrows indicate optimal values of 
w, and we report all choices of w in Table [T] 
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FIG. 4. (color online) The pair correlation function g(r) for 32 
deuterium atoms at T — 5000 K, computed using PIMC with 
augmented nodes. The contribution to g(r) from an isolated 
molecule at 300 K is shown for comparison, illustrating the 
compression of the molecules. 



partition function decreases for larger values of w, indi- 
cating an optimal mixture of atomic Is and free particle 
density matrices occurs near w = 5. At T — 10 000 
a value of w around 3 or larger increases the partition 
function by more than a factor of 20. The partition func- 
tion is relatively insensitive to w for all w > 3, indicating 
that the Is orbitals dominate the nodal surface. 

Such a systematic approach for optimizing a given 
nodal model from first principles has not been possible 
before within the reference slice approach, due to the fact 
that the nodes are time-dependent and require a com- 
plete set of terms over the full range of inter-slice tem- 
peratures from MT down to 2T. We note that this is 
analogous to optimizing ground nodes directly diffusion 
Monte Carlo energies, which is the most rigorous ground 
state fixed-node strategy [22J, [23] . 

To allow direct comparison with traditional reference 
slice PIMC, we performed calculations of the hydrogen 
Hugoniot using the same densities as Ref. |7J We ran 
periodic simulation with N = 32 and N = 64 deu- 
terium atoms and for both densities with r s = 2.00 and 
r s = 1.86. Following Ref. [7J we identify the Hugoniot 
curve by linear interpolation. The two densities accu- 



rately bracket the Hugoniot curve at higher pressures, 
but was found to introduce an extrapolation error at the 
lowest pressures, making the statistical error bars shown 
in Fig. [T] a lower bound only. Following the guidance in 
Ref. [7J on the importance of small time steps, we used 
Ar/h = 0.02 Ha -1 , a time step at least twice as small as 
Ref. [7) Our results are summarized in Fig. [T] and Table [I] 

In the high temperature and pressure limit, our Hugo- 
niot approaches the correct free particle limit of four- 
fold compression. We see excellent agreement with 
past PIMC results at high temperatures and pressures 
[7J, where permutations are infrequent and free particle 
nodes accurately describe a dense strongly interacting 
plasma. An augmented node description starts becoming 
necessary below 62 500 K, resulting in lower free energy 
than free particle nodes alone. Using augmented nodes, 
we see the Hugoniot curving towards lower compression 
in the 1 Mbar regime, before eventually coming closer 
to gas-gun and DFT low pressure results [31 [23]. We 
note that the antipodal slice approximation also speeds 
up simulations considerably, as we find that parallel sim- 
ulations run more than four times faster in direct com- 
parisons to reference-slice fixed-node simulations. 

We note that these results differ somewhat from the 
lowest pressure results of Ref. [7], which used a varia- 
tional approach to optimize the parameters in a density 
matrix assumed to be a Slater determinant of Gaussians. 
The variational nodes are shown to give lower energies 
than the free particle density matrix in Ref. [7J however, 
no formal proof is given that variational nodes used there 
result in the lowest free energy. Further investigation is 
needed to resolve this discrepancy. We note that using 
Is orbitals in our augmented node results is also a crude 
model that can be enhanced by adding higher orbitals. 

Our results do not seem to suffer from finite size effects 
at higher temperatures. This is not the case at low tem- 
peratures. We comment on the exceptional discrepancy 
at T= 15 625 K. This is where DFT-MD shows some dif- 
ficulty probing the simultaneous molecular bonding and 
dissociations that take place in this range. DFT which 
can be notorious in describing bond lengths and band gap 
closure could be exaggerating a high compression trend. 
Long distance correlations and delocalized electron paths 
are observed in our PIMC simulations. This may suggest 
the need for larger PIMC simulation box size and more 
atoms to account for these correlations. 

We note that in the direct PIMC approach by [IB] , the 
PIMC simulations did not seem to converge to an equilib- 
rium state around 10 000 K, and a plasma transition into 
partially ionized deuterium was assumed to take place in 
that region. A notable configurational expression of this 
is the nucleii form metallic clusters (droplets) with highly 
delocalized electrons over them [THJ. We do not observe 
this behavior in our simulation. 

Our data are within error bars of the magnetically 
driven flyer [15j and convergent explosive |16j experi- 
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TABLE I. PIMC calculation of deuterium Hugoniot and compressibilities with respect to an initial cryogenic state of 0.171 g/cc 
at T = 19.6 K. A small time step of r _1 = 1.536 x 10 7 K (Ar/fi = 0.02 Ha~\ twice smaller than Ref. 0) is used. 
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ments. If the recent correction on the equation of state of 
the quartz standard 25] is taken into account, the laser 
data shifts to lower compressions indicating that Hydro- 
gen is not as compressible as the older data suggested. 

In summary, we present a new time-independent for- 
malism for fixed-node fermionic PIMC which is signif- 
icantly more efficient and free of the ergodicity issues 
that plague the traditional reference slice approach. A 
clear strategy for creating nodal models with lower free 
energy emerges from this approach, which we anticipate 
will open the door to systems previously inaccessible to 
PIMC, including elements with Z > 2. We present 
updated results on the hydrogen Hugoniot which show 
a clear trend toward the lower compression and closer 
agreement with Kerley's equation of state [25] . 
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